1. Импорт данных и первичный анализ

data <- readRDS("life_expectancy_data.RDS")

Описание данных

data_info <- describe(data, na.rm = TRUE, skew = FALSE, ranges = TRUE)
data_info
##                                         vars   n         mean           sd
## Country*                                   1 195 9.800000e+01 5.644000e+01
## Year                                       2 195 2.019000e+03 0.000000e+00
## Gender*                                    3 195 1.000000e+00 0.000000e+00
## Life expectancy                            4 195 7.552000e+01 7.690000e+00
## Unemployment                               5 195 8.600000e+00 6.970000e+00
## Infant Mortality                           6 195 1.961000e+01 1.693000e+01
## GDP                                        7 195 4.659960e+11 1.933272e+12
## GNI                                        8 195 4.864381e+11 1.956716e+12
## Clean fuels and cooking technologies       9 195 6.598000e+01 3.632000e+01
## Per Capita                                10 195 1.682098e+04 2.444060e+04
## Mortality caused by road traffic injury   11 195 1.706000e+01 1.037000e+01
## Tuberculosis Incidence                    12 195 1.037500e+02 1.342800e+02
## DPT Immunization                          13 195 8.799000e+01 1.240000e+01
## HepB3 Immunization                        14 195 8.676000e+01 1.273000e+01
## Measles Immunization                      15 195 8.731000e+01 1.319000e+01
## Hospital beds                             16 195 3.000000e+00 2.370000e+00
## Basic sanitation services                 17 195 7.738000e+01 2.777000e+01
## Tuberculosis treatment                    18 195 7.757000e+01 1.705000e+01
## Urban population                          19 195 5.912000e+01 2.329000e+01
## Rural population                          20 195 4.088000e+01 2.329000e+01
## Non-communicable Mortality                21 195 1.705000e+01 7.080000e+00
## Sucide Rate                               22 195 4.800000e+00 3.780000e+00
## continent*                                23 195 2.670000e+00 1.310000e+00
##                                                  min          max        range
## Country*                                        1.00 1.950000e+02 1.940000e+02
## Year                                         2019.00 2.019000e+03 0.000000e+00
## Gender*                                         1.00 1.000000e+00 0.000000e+00
## Life expectancy                                55.49 8.810000e+01 3.260000e+01
## Unemployment                                    0.18 3.644000e+01 3.626000e+01
## Infant Mortality                                1.40 7.580000e+01 7.440000e+01
## GDP                                     188391770.64 2.143322e+13 2.143304e+13
## GNI                                     375392118.22 2.170865e+13 2.170827e+13
## Clean fuels and cooking technologies            0.00 1.000000e+02 1.000000e+02
## Per Capita                                    228.21 1.758139e+05 1.755857e+05
## Mortality caused by road traffic injury         0.00 6.460000e+01 6.460000e+01
## Tuberculosis Incidence                          0.00 6.540000e+02 6.540000e+02
## DPT Immunization                               35.00 9.900000e+01 6.400000e+01
## HepB3 Immunization                             35.00 9.900000e+01 6.400000e+01
## Measles Immunization                           37.00 9.900000e+01 6.200000e+01
## Hospital beds                                   0.20 1.371000e+01 1.351000e+01
## Basic sanitation services                       8.63 1.000000e+02 9.137000e+01
## Tuberculosis treatment                          0.00 1.000000e+02 1.000000e+02
## Urban population                               13.25 1.000000e+02 8.675000e+01
## Rural population                                0.00 8.675000e+01 8.675000e+01
## Non-communicable Mortality                      4.40 4.370000e+01 3.930000e+01
## Sucide Rate                                     0.30 3.010000e+01 2.980000e+01
## continent*                                      1.00 5.000000e+00 4.000000e+00
##                                                   se
## Country*                                4.040000e+00
## Year                                    0.000000e+00
## Gender*                                 0.000000e+00
## Life expectancy                         5.500000e-01
## Unemployment                            5.000000e-01
## Infant Mortality                        1.210000e+00
## GDP                                     1.384445e+11
## GNI                                     1.401234e+11
## Clean fuels and cooking technologies    2.600000e+00
## Per Capita                              1.750230e+03
## Mortality caused by road traffic injury 7.400000e-01
## Tuberculosis Incidence                  9.620000e+00
## DPT Immunization                        8.900000e-01
## HepB3 Immunization                      9.100000e-01
## Measles Immunization                    9.400000e-01
## Hospital beds                           1.700000e-01
## Basic sanitation services               1.990000e+00
## Tuberculosis treatment                  1.220000e+00
## Urban population                        1.670000e+00
## Rural population                        1.670000e+00
## Non-communicable Mortality              5.100000e-01
## Sucide Rate                             2.700000e-01
## continent*                              9.000000e-02
na_count_cols <- colSums(is.na(data))
na_count_cols #Сколько N/A в каждом столбце
##                                 Country                                    Year 
##                                       0                                       0 
##                                  Gender                         Life expectancy 
##                                       0                                       0 
##                            Unemployment                        Infant Mortality 
##                                       0                                       0 
##                                     GDP                                     GNI 
##                                       0                                       0 
##    Clean fuels and cooking technologies                              Per Capita 
##                                       0                                       0 
## Mortality caused by road traffic injury                  Tuberculosis Incidence 
##                                       0                                       0 
##                        DPT Immunization                      HepB3 Immunization 
##                                       0                                       0 
##                    Measles Immunization                           Hospital beds 
##                                       0                                       0 
##               Basic sanitation services                  Tuberculosis treatment 
##                                       0                                       0 
##                        Urban population                        Rural population 
##                                       0                                       0 
##              Non-communicable Mortality                             Sucide Rate 
##                                       0                                       0 
##                               continent 
##                                       0

 

 

 

2. Интерактивный plotly график для двух нумерических переменных (), раскрашенный по колонке континента, на котором расположена страна

plot_ly(
  data[data$`Life expectancy` != 0,],
  x = ~ `Life expectancy`,
  color = ~continent,
  type = "box"
)
plot_ly(
  data[data$`Sucide Rate` != 0,],
  x = ~ `Sucide Rate`,
  color = ~continent,
  type = "box"
)

 

 

 

3. Cравнение распределений колонки Life expectancy между группами стран Африки и Америки. Вид статистического теста определите самостоятельно. Визуализируйте результат через библиотеку rstatix.

data_filtered <- data %>%
  filter(continent %in% c("Africa", "Americas"))

test_data_filtered <- data_filtered %>%
  t_test(`Life expectancy` ~ continent,
         data = .,
         paired = FALSE,
         var.equal = TRUE)

print(test_data_filtered)
## # A tibble: 1 × 8
##   .y.             group1 group2      n1    n2 statistic    df        p
## * <chr>           <chr>  <chr>    <int> <int>     <dbl> <dbl>    <dbl>
## 1 Life expectancy Africa Americas    52    38     -11.4    88 4.12e-19

Поскольку p-value ниже заданного значения 0.05, то мы отвергаем Н0 и можем утверждать, что продолжительность жизни отличается в странах Африки и Америки.

sum_stats <- data_filtered %>%
  get_summary_stats(`Life expectancy`, type = "mean_sd")
sum_stats
## # A tibble: 1 × 4
##   variable            n  mean    sd
##   <fct>           <dbl> <dbl> <dbl>
## 1 Life expectancy    90  71.4  8.05
library(ggpubr)
ggqqplot(data_filtered[data_filtered$`Life expectancy`, ], 
         x = "Life expectancy", facet.by = "continent")

#Данные не выглядят нормально распределенными, используем тест Манна-Уитни вместо Т-теста
stat.test <- data_filtered[data_filtered$`Life expectancy`, ] %>% 
  wilcox_test(`Life expectancy` ~ continent) %>% 
  add_xy_position(x = "Life expectancy")
stat.test
## # A tibble: 1 × 11
##   .y.       group1 group2    n1    n2 statistic        p y.position groups  xmin
##   <chr>     <chr>  <chr>  <int> <int>     <dbl>    <dbl>      <dbl> <name> <dbl>
## 1 Life exp… Africa Ameri…    47    43        72 3.25e-14       84.1 <chr>     NA
## # ℹ 1 more variable: xmax <dbl>
library(rstatix)

ggboxplot(
  data_filtered[data_filtered$`Life expectancy`, ], 
  x = "continent", y = "Life expectancy", 
  ylab = "Life expectancy", xlab = "Continent", 
  add = "jitter"
  ) + 
  labs(subtitle = get_test_label(stat.test, detailed = FALSE)) + 
  stat_pvalue_manual(stat.test, tip.length = 0) 

 

 

 

4. Сделайте новый датафрейм, в котором оставите все численные колонки кроме Year. Сделайте корреляционный анализ этих данных. Постройте два любых типа графиков для визуализации корреляций.

data_filtered_2 <- data %>%
  select_if(is.numeric) %>%
  select(-Year, -`Life expectancy`)

data_f2_cor <- cor(data_filtered_2) 
library(corrplot)
## corrplot 0.92 loaded
corrplot(data_f2_cor, method = 'number', number.cex = 0.8, tl.cex = 0.7)

corrplot(data_f2_cor, method = "color", type = "lower", 
         addCoef.col = "grey29", diag = FALSE,
         cl.pos = "r", tl.col = "grey11",
         col = COL2('RdBu', 10),
         tl.cex = 0.8)

 

 

 

5. Постройте иерархическую кластеризацию на этом датафрейме.

library(factoextra)
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
#Стандартизация данных
data_scaled <- scale(data_filtered_2)

#Матрица дистанций
distances_data_scaled <- dist(data_scaled, 
                        method = "euclidean"
                        )
as.matrix(distances_data_scaled)[1:6, 1:6]
##          1        2        3        4        5        6
## 1 0.000000 7.390287 6.144642 4.404280 6.468405 7.723136
## 2 7.390287 0.000000 2.610829 7.637350 3.346255 3.630919
## 3 6.144642 2.610829 0.000000 6.049814 4.350330 3.456356
## 4 4.404280 7.637350 6.049814 0.000000 7.886278 6.853670
## 5 6.468405 3.346255 4.350330 7.886278 0.000000 4.960145
## 6 7.723136 3.630919 3.456356 6.853670 4.960145 0.000000
#Дендрограмма кластеров
data_scaled_hc <- hclust(d = distances_data_scaled, 
                        method = "ward.D2")

#Визуализация
fviz_dend(data_scaled_hc, 
          cex = 0.1) # cex() - размер лейблов

 

 

 

6. Сделайте одновременный график heatmap и иерархической кластеризации. Содержательно интерпретируйте результат.

library(pheatmap)


plot3 <- pheatmap(data_scaled, 
                  show_rownames = FALSE, 
                  clustering_distance_rows = distances_data_scaled,
                  clustering_method = "ward.D2", 
                  cutree_rows = 5,
                  cutree_cols = length(colnames(data_scaled)),
                  angle_col = 45, 
                  main = "Heatmap + hierarchical clustering")

plot3

В данном графике дендрограмма сверху показывает иерархическую кластеризацию переменных. Кластеры отображают какие переменные схожи, а высота соединения показывает степень различия между группами переменных. Два больших кластера. Первый с переменными, связанными с уровнем безработицы, туберкулезом и его лечением, разными Mortality. Во втором кластере переменные GDP и GNI являются похожими и стоят отдельно от всех остальных переменных в этом кластере. Переменные связанные с иммунизацией сгруппировались в один схожий кластер. Вертикальная дендрограмма (сбоку) отражает иерархическую кластеризацию наблюдений (объектов). Здесь выделяется 5 групп (одна очень маленькая с высокими стандартизированными значениями переменных GDP, GNI). Эти 5 групп, вероятно, более-менее отражают группы по континентам.

counts <- table(data$continent)

# Вывод результатов
print(counts)
## 
##   Africa Americas     Asia   Europe  Oceania 
##       52       38       42       48       15

 

 

 

7. Проведите PCA анализ на этих данных. Проинтерпретируйте результат.

library(FactoMineR)

data_full.pca <- prcomp(data_filtered_2, 
                        scale = T) 
summary(data_full.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.6027 1.4816 1.3933 1.16723 1.08021 0.96332 0.92844
## Proportion of Variance 0.3764 0.1220 0.1078 0.07569 0.06483 0.05155 0.04789
## Cumulative Proportion  0.3764 0.4983 0.6061 0.68184 0.74666 0.79822 0.84611
##                            PC8     PC9    PC10    PC11    PC12    PC13   PC14
## Standard deviation     0.84336 0.69074 0.68425 0.58700 0.54757 0.42231 0.3550
## Proportion of Variance 0.03951 0.02651 0.02601 0.01914 0.01666 0.00991 0.0070
## Cumulative Proportion  0.88562 0.91213 0.93814 0.95728 0.97394 0.98385 0.9909
##                           PC15    PC16    PC17      PC18
## Standard deviation     0.34451 0.20224 0.07146 3.138e-16
## Proportion of Variance 0.00659 0.00227 0.00028 0.000e+00
## Cumulative Proportion  0.99744 0.99972 1.00000 1.000e+00
fviz_eig(data_full.pca, addlabels = T, ylim = c(0, 40))

 

Первые две компоненты объясняют около 50% дисперсии.

fviz_pca_var(data_full.pca, col.var = "contrib", labelsize = 3)

fviz_pca_var(data_full.pca, 
             select.var = list(contrib = 5), # Задаём число здесь 
             col.var = "contrib")

fviz_contrib(data_full.pca, choice = "var", axes = 1, top = 24) # 1

fviz_contrib(data_full.pca, choice = "var", axes = 2, top = 24) # 2

fviz_contrib(data_full.pca, choice = "var", axes = 3, top = 24) # 3

 

Infant Mortality and Clean fuels and cooking technologies стали самыми важными переменными с точки зрениях их вариации в PC1. Иммунизации стали тремя самыми важными переменными с точки зрениях их вариации в PC2.

 

 

 

8.Постройте biplot график для PCA. Раскрасьте его по значениям континентов. Переведите его в plotly. Желательно, чтобы при наведении на точку, вы могли видеть название страны.

ggbiplot(data_full.pca, 
         scale=0, alpha = 0.1) + 
  theme_minimal()

Более осмысленным biplot становится при использовании кластерных методов, с помощью которых мы можем разделить наблюдения на группы. Посмотрим, наблюдается ли разница между группами по continents:

# Сделаем корректные данные для группировки по continents.

data_filtered_2_with_continent <- data_filtered_2 %>%
  mutate(continent = data$continent)



# Визуализируем с группировкой по diabetes (для этого переменную нужно сделать фактором)
ggbiplot(data_full.pca, 
         scale=0, 
         groups = as.factor(data_filtered_2_with_continent$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal()

pca_df <- as.data.frame(data_full.pca$x)
pca_df$country <- data$Country
pca_df$continent <- as.factor(data_filtered_2_with_continent$continent)

p <- plot_ly(pca_df, x = ~PC1, y = ~PC2, text = ~country, color = ~continent, 
             type = 'scatter', mode = 'markers',
             marker = list(opacity = 0.7)) %>%
  layout(title = 'PCA Plot',
         xaxis = list(title = 'PC1'),
         yaxis = list(title = 'PC2'))

p

 

 

 

9.Дайте содержательную интерпретацию PCA анализу.

Африка сильно отличается от Европы, Америки и Азии. Океания самая самая “разбросанная”, но в ней и меньше всего объектов. Infant Mortality and Clean fuels and cooking technologies стали самыми важными переменными с точки зрениях их вариации в PC1, а три вида иммунизации стали тремя самыми важными переменными с точки зрениях их вариации в PC2.  

 

 

10. Сравните результаты отображения точек между алгоритмами PCA и UMAP.

Tidymodels approach

#install.packages("embed")
library(tidymodels)
library(embed)

umap_prep <- recipe(~., data = data_filtered_2) %>% # "техническая" строка, нужная для работы фреймворка tidymodels
  step_normalize(all_predictors()) %>% # нормируем все колонки
  step_umap(all_predictors()) %>%  # проводим в UMAP. Используем стандартные настройки. Чтобы менять ключевой параметр (neighbors), нужно больше погружаться в машинное обучение
  prep() %>%  # "техническая" строка, нужная для работы фреймворка tidymodels. Мы выполняем все степы выше 
  juice() # Финальная строка - приводим результаты UMAP к стандартизированному датасету

Визуализиуем два первых измерения UMAP и добавим информацию о возрастных группах и диабет-статусе:

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + #  # можно добавить раскраску 
  geom_point(aes(color = as.character(data_filtered_2_with_continent$continent)), 
             alpha = 0.7, size = 2) +
  labs(color = NULL) 

 

Отображение точек в целом схоже в PCA и в UMAP.

 

 

 

11. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Удалите 5 случайных колонок. Проведите PCA анализ. Повторите результат 3 раза. Наблюдаете ли вы изменения в куммулятивном проценте объяснённой вариации? В итоговом представлении данных на биплотах? С чем связаны изменения между тремя PCA?

set.seed(123)
columns_to_remove_1 <- sample(names(data_filtered_2), 5)
data_filtered_3 <- select(data_filtered_2, -all_of(columns_to_remove_1))
columns_to_remove_1
## [1] "Urban population"       "Tuberculosis treatment" "GDP"                   
## [4] "HepB3 Immunization"     "Infant Mortality"
columns_to_remove_2 <- sample(names(data_filtered_2), 5)
data_filtered_4 <- select(data_filtered_2, -all_of(columns_to_remove_2))
columns_to_remove_2
## [1] "Measles Immunization"                
## [2] "Clean fuels and cooking technologies"
## [3] "GNI"                                 
## [4] "Tuberculosis treatment"              
## [5] "Per Capita"
columns_to_remove_3 <- sample(names(data_filtered_2), 5)
data_filtered_5 <- select(data_filtered_2, -all_of(columns_to_remove_3))
columns_to_remove_3
## [1] "Clean fuels and cooking technologies"
## [2] "DPT Immunization"                    
## [3] "Basic sanitation services"           
## [4] "GDP"                                 
## [5] "Tuberculosis Incidence"
data_full.pca_3 <- prcomp(data_filtered_3, 
                        scale = T) 
fviz_eig(data_full.pca_3, addlabels = T, ylim = c(0, 40))

data_full.pca_4 <- prcomp(data_filtered_4, 
                        scale = T) 
fviz_eig(data_full.pca_4, addlabels = T, ylim = c(0, 40))

data_full.pca_5 <- prcomp(data_filtered_5, 
                        scale = T) 
fviz_eig(data_full.pca_5, addlabels = T, ylim = c(0, 40))

 

Кумулятивный процент объясненной вариации для PC1 и PC2 при выбрасывании трех случайных переменных составил 50.8%, 54.2% и 47.6% (по сравнению с 49.8% при предыдущем анализе)

library(cowplot)

pl3 <- ggbiplot(data_full.pca_3, scale=0, alpha = 0.1) + theme_minimal()
pl4 <- ggbiplot(data_full.pca_4, scale=0, alpha = 0.1) + theme_minimal()
pl5 <- ggbiplot(data_full.pca_5, scale=0, alpha = 0.1) + theme_minimal()
plot_grid(pl3, pl4, pl5, ncol = 3)

 

При выбрасывании 5 переменных достаточно сильно изменилось отображение биплотов, хотя и некоторые переменные остались примерно в том же направлении (например, Hospital beds).

 

 

 

12. Давайте самостоятельно увидим, что снижение размерности – это группа методов, славящаяся своей неустойчивостью. Создайте две дамми-колонки о том: (1) принадлежит ли страна к африканскому континенту, (2) Океании. Проведите PCA вместе с ними, постройте биплоты. Проинтерпрейтируйте результат. Объясните, почему добавление дамми-колонок не совсем корректно в случае PCAнашего типа.

data_filtered_6 <- data_filtered_2_with_continent %>%
  mutate(`is Africa` = as.integer(continent == "Africa"),
         `is Oceanina` = as.integer(continent == "Oceania")) %>%
  select(-continent)
data_full.pca_6 <- prcomp(data_filtered_6, 
                        scale = T) 
fviz_eig(data_full.pca_6, addlabels = T, ylim = c(0, 40))

47.2% вариативности объясняют PC1 и PC2.

pl6 <- ggbiplot(data_full.pca_6, scale=0, alpha = 0.1) + theme_minimal()
pl6

# Сделаем корректные данные для группировки по continents.

data_filtered_6_with_continent <- data_filtered_6 %>%
  mutate(continent = data$continent)



# Визуализируем с группировкой по diabetes (для этого переменную нужно сделать фактором)
ggbiplot(data_full.pca_6, 
         scale=0, 
         groups = as.factor(data_filtered_6_with_continent$continent), 
         ellipse = T,
         alpha = 0.2) +
  theme_minimal()

fviz_pca_var(data_full.pca_6, col.var = "contrib", labelsize = 3)

fviz_pca_var(data_full.pca_6, 
             select.var = list(contrib = 5), # Задаём число здесь 
             col.var = "contrib")

fviz_contrib(data_full.pca_6, choice = "var", axes = 1, top = 24) # 1

fviz_contrib(data_full.pca_6, choice = "var", axes = 2, top = 24) # 2

fviz_contrib(data_full.pca_6, choice = "var", axes = 3, top = 24) # 3

 

В целом, результаты не сильно изменились в сравнении с результатами без дамми-колонок, однако привели к изменению порядка состава главных компонент. В Dim-1 новая дамми-переменная ‘is Africa’ стала 6-й по вкладу в вариативность (около 7%), а новая дамми-переменная ‘is Oceania’ сильно не повлияла (во все три Dim ее влияние около 0%), что связано с большим числом (большей вариативности) объектов из Африки. Также, если вариативность дамми-переменных была бы значительна по сравнению со всеми остальными и было бы необходимо включить такие дамми-переменные в анализ, можно было бы стандартизировать все остальные переменные, чтобы дамми-переменные не искажали масштаб. В нашем случае целью было исследовать взаимосвязи между количественными переменными при помощи евклидовых расстояний, поэтому включение переменных ‘is Africa’ и ‘is Oceania’, которые по факту являются категориальными не лучшее решение.